import matplotlib.pyplot as plt
import numpy as np
import rasterio
import os

nodata = None
mean = list()
path = './mean_npp/'
for file in sorted(os.listdir(path)):
    with rasterio.open(os.path.join(path,file)) as src:
        nodata = src.nodata
        print("文件路径:", src.name)
        print("数据类型:", src.dtypes[0])  # 查看第一个波段的数据类型
        print("波段数:", src.count)  # 波段数量
        print("影像形状 (行, 列):", src.shape)  # (高度, 宽度)
        print("分辨率:", src.res)  # 像元大小 (X 方向, Y 方向)
        print("投影坐标系 (CRS):", src.crs)  # 坐标参考系统
        print("仿射变换 (GeoTransform):", src.transform)  # 地理仿射变换矩阵
        print("NoData 值:", src.nodata)  # 无数据值
        print("元数据:", src.meta)  # 影像的所有元数据
        print("Scale Factor:", src.scales)  # 检查比例因子
        print("Offset:", src.offsets)  # 检查偏移量
        image = src.read(1)  # 读取第一个波段
        print(np.unique(image))
        image_except_nodata = image[image!=nodata]#*src.scales
        mean.append(image_except_nodata.mean())
# 统计 NoData 之外的值

mean_npp = np.array(mean)
print('mean:{}'.format(mean_npp))
years = np.arange(2000, 2000 + len(mean_npp))

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# 上图：年NPP值和趋势线
ax1.plot(years, mean_npp, marker='s', label='年NPP值')
z = np.polyfit(years, mean_npp, 1)
p = np.poly1d(z)
ax1.plot(years, p(years), "r-", label='趋势线')
ax1.set_ylabel('NPP变化 (gc·m⁻²·a⁻¹)')
ax1.legend()
ax1.text(0.05, 0.95, f'y={z[0]:.4f}x+{z[1]:.2f}\n$R^2$={np.corrcoef(years, mean_npp)[0, 1]**2:.2f}', 
         transform=ax1.transAxes, fontsize=12, verticalalignment='top')

# 下图：NPP偏离值
deviations = mean_npp - p(years)
ax2.bar(years, deviations, label='NPP偏离值')
ax2.set_ylabel('Annual variations of NPP (gc·m⁻²·a⁻¹)')
ax2.set_xlabel('年份')
ax2.legend()

plt.tight_layout()
plt.show()
# print('平均值:{}'.format(mean))

